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ABSTRACT 


In  combination  with  an  appropriate  background  ocean  model,  simulations  of  internal 
waves  are  necessary  to  create  realistic  scenarios  for  use  in  acoustic  propagation  studies. 
Our  approach  employs  contemporary  understanding  of  internal  wave  statistics  with  a  sim¬ 
ulator  that  produces  two-dimensional  (range/depth)  slices  of  internal  wave  displacements. 
The  method  produces  statistical  realizations  for  acoustic  studies  that  are  true  to  ocean 
measurements  with  nominal  computational  cost. 
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1  Introduction 


This  report  describes  a  numerical  simulator  that  produces  realizations  of  internal  wave  dis¬ 
placements  for  use  in  studies  of  long-range  acoustic  propagation.  The  goals  of  any  simulator 
are  twofold:  first  to  produce  simulations  that  are  realistic,  and  second,  to  do  so  efficiently. 

There  are  a  variety  of  numerical  simulators  for  internal  waves.  Like  the  one  described 
here,  most  are  spectral  models.  These  models  commonly  utilize  the  formulation  for  the 
Garrett  and  Munk  (GM)  spectrum  and  in  particular,  the  formulation  based  on  a  simplified 
Wentzel-Kramers-Brillouin  (WKB)  approximation  [6].  In  shallow  water  or  near  the  surface, 
these  WKB  approximations  are  not  valid.  Taking  a  non- WKB  approach  improves  simulation 
realism  with  little  cost,  and  provides  an  algorithm  valid  in  both  deep  and  shallow  water. 

For  several  years,  the  simulator  described  as  a  starting  field  in  Winters  and  D’Asaro 
[7]  (hereinafter  WD97)  was  employed  to  model  situations  involving  relatively  short-range 
acoustic  transmission  experiments,  i.e.,  less  than  50  km.  The  fields  may  be  incremented 
in  time  consistent  with  the  internal  wave  dispersion  relationship  and  the  method  is  almost 
free  of  the  WKB  limitation.  For  shallow  water  and  short-range  situations,  their  approach 
was  extended  to  relax  this  limitation.  However,  because  the  code  is  three-dimensional,  the 
approach  is  computationally  intensive.  And,  to  minimize  memory  usage,  strict  limits  are 
imposed  on  the  smallest  and  largest  wave  numbers  that  may  be  represented. 

An  attempt  was  made  to  use  3-D  fields  from  the  WD97  technique  for  long-range  prop¬ 
agation  studies.  Two-dimensional  depth-range  fields  were  obtained  by  slicing  through  the 
3-D  field  along  an  angle,  then  interpolating  the  field  onto  the  desired  grid.  This  allows 
continuous  slices  to  any  range  to  be  formed.  In  WD97,  2-D  horizontal  fast  Fourier  trans¬ 
forms  (FFTs)  are  taken  at  each  depth  so  the  fields  are  periodic  about  the  boundaries. 
Unfortunately,  periodicities  related  to  the  size  of  the  box  inevitably  are  produced  in  the 
output  fields.  So  WD97  produces  reasonably  realistic  fields  within  its  domain.  However,  it 
is  computationally  intensive  and  when  extended  to  long  range,  the  resulting  fields  are  not 
realistic.  The  WD97  approach  was  not  designed  for  long-range  propagation  studies. 

It  was  realized  that  by  sampling  in  kx  and  k  =  \J~k^~+~k^  rather  than  the  traditional 
kx  and  ky,  only  a  moderate  number  of  k- values  are  required.  But  a  trade-off  is  made:  the 
problem  reduces  to  a  single  Fourier  transform  in  x  providing  a  field  at  one  (or  a  few)  y- 
values.  The  resulting  internal  wave  simulator  is  very  efficient  and  meets  our  requirements 
of  producing  displacement  fields  that  (1)  have  nreganreter  or  greater  periodicities,  (2)  are 
not  constrained  by  the  WKB  approximation,  and  (3)  contain  a  broad  range  of  vertical 
and  horizontal  scales.  Hence  the  simulator  generates  realistic  fields  to  very  long  ranges 
(megameters).  In  section  II,  the  formulation  is  described.  In  section  III,  examples  and 
results  from  tests  are  shown.  An  appendix  presents  a  pseudo-code  for  reference  by  those 
wishing  copies  of  the  software. 
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2  Method 


Much  of  our  formalism  develops  from  Henyey  et  al.  [4]  (hereinafter  H97).  The  numerical 
simulator  follows  the  presentation  in  WD97,  where  a  general  method  for  a  five-variable 
simulation  is  presented.  Here,  obtaining  sound  speed  variations  requires  only  the  internal 
wave  displacement  field.  Depict  this  field  in  two  dimensions  as  £(x,  y,  Z',t).  The  third 
dimension  y  =  0.  Note  that  throughout  this  discussion,  time  is  considered  a  parameter. 
For  convenience  consider  the  displacement  field  in  terms  of  the  component  horizontal  wave 
number  kx  ,  where  x  and  kx  are  Fourier  transform  pairs,  then 

r  M 

C (kx,  y,z)  =  /  eiyky  J2{G-e-^kt  +  G+e^jcpj^dky  (1) 

J  3  =  1 


which  is  analogous  to  Eq.  18  in  WD97. 

We  evaluate  this  expression  at  a  given  y  =  y$.  An  inverse  Fourier  transform  at  each 
depth  provides  ((x,  yo,  z;  t).  In  an  ocean  modeled  as  horizontally  infinite  in  both  directions, 
G-  and  G+  are  generalized  functions  of  a  continuous  wave  vector.  From  this  point  we  will, 
as  an  approximation,  treat  them  as  ordinary  functions  of  a  discrete  wave  vector.  Each  of  the 
discrete  wave  vectors  is  representative  of  a  small  area  of  continuous  wave  vector  surrounding 
the  discrete  value.  These  areas  fill  up  the  whole  space  without  overlapping.  The  value  of  G_ 
or  G+  at  one  of  the  discrete  points  is  the  average  of  the  corresponding  generalized  function 
over  the  area  surrounding  that  point. 


For  a  finite  depth  ocean,  there  is  a  decomposition  in  the  wave  field  in  vertical  modes 
< Pj:k(z ),  where  mode  number  j  =  1,2,  The  modes  y>j,k(z )  and  their  frequencies 

are  functions  of  the  horizontal  wave  number  k,  where  k2  =  k2  +  k2  (cf.  WD97  Eq.  8 
ff.).  The  modes  and  frequencies  Ujk  are  obtained  by  a  method  that  has  been  passed  from 
researcher  to  researcher  for  many  years.  We  map  the  problem  onto  finding  the  eigenvalues 
and  eigenvectors  of  a  symmetric,  tri-diagonal  matrix.  In  terms  of  eigenfunctions  = 
y/ n2,  —  f2tp  ,  Eq.  8  of  WD97  for  the  vertical  modes  may  be  written  as 


Vn2  ~  f 


-;(k2~d2z) 


\/n2  -  f 


=  \lf) 


(2) 


The  eigenvalues 


X  =  k2/(co2-f2) 


(3) 


and  eigenfunctions  are  obtained  rising  a  tri-diagonal  mode  solver, 
mode  normalization 

ro 


' —Lz 


{n2  ~  f2)Vj,k{z)Vj',k(z)  dz  =  5. 


n 


In  this  mapping,  the 


(4) 


and  boundary  conditions,  Lz)  =  0)  =  0  are  preserved. 
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To  obtain  the  displacement  field  in  two  dimensions  for  a  given  depth  as  a  function  of 
horizontal  wave  number,  there  is  a  summation  over  j  and  over  k.  This  differs  from  WD97. 
The  reduction  to  two  dimensions  requires  us  to  define  ky  in  terms  of  k  and  kx. 

The  wave  amplitudes  are  determined  by  G-  and  G. |_.  These  are  complex  random  de¬ 
viates  generated  separately  to  allow  independent  right-  and  left-going  waves.  G-  and  £?+ 
have  zero  mean,  with  a  variance  given  by  the  internal  wave  spectrum  defined  below  and 
by  the  wave  number  area  represented  by  one  wave  vector.  The  variance  is  proportional  to 
1/area  because  of  the  averaging  of  delta-correlated  quantities. 

The  technique  defined  in  Eq.  1  is  common.  Realizations  are  obtained  via  Monte  Carlo 
method  on  the  weights,  G-  and  G+.  A  displacement  realization  is  generated  in  the  wave 
number  domain,  with  amplitude  determined  by  the  model  spectrum.  Note  that  the  depth 
dependence  is  specified  by  the  amplitudes  of  the  vertical  modes.  (This  property  is  one 
step  in  freeing  us  from  the  WKB  approximation.)  At  this  stage,  WD97  uses  a  2-D  inverse 
Fourier  transform  at  each  depth  and  obtains  the  field  in  the  3-D  spatial  domain. 

Our  method  uses  the  full  spectral  model.  This  is  reduced  to  one  dimension  in  the  hori¬ 
zontal  by  summing  over  circles  or  annuli  in  horizontal  wave  number  space.  Once  generated 
in  the  wave  number  domain,  only  a  one- dimensional  inverse  FFT  is  needed  to  obtain  the 
field  in  range.  The  result  is  an  efficient  and  accurate  algorithm. 


Consider  the  differential  displacement  variance  in  frequency,  mode  number,  and  hori¬ 
zontal  wave  number 

S(u,j)dudj ^  (5) 

Though  the  mode  number  increment,  dj,  equals  1,  we  leave  it  in  for  completeness.  We 
assume  the  spectrum  and  the  vertical  modes  are  independent  of  horizontal  direction  (i.e., 
the  field  is  horizontally  isotropic).  In  Eq.  5,  0  is  the  direction  of  the  horizontal  wave  number. 
Writing  the  variance  in  terms  of  increments  in  horizontal  wave  numbers  gives 

S(u,j)du;dj ^  =  S^kdk' '  dddkxdky  (6) 

Recognizing  the  group  speed,  vg  =  duo/dk,  and  folding  the  ky  integral  to  the  positive  ky 
side,  the  differential  variance  becomes 


var  = 

In  terms  of  the  annulus  ki  <  k  <  kf, 


S{u,j)vg 
TT  k 


djdkxd\ky 


dk  =  kf  —  kj. 

d\ky\  =  yjkj-k^-yj kf  -  k2x 


(7) 

(8) 

(9) 


At  this  point,  our  derivation  has  diverged  from  H97.  Finally,  using  Eq.  6,  with 
an  expected  value, 


(|G+|2)  +  (|G_|2)  =  S{uJ^k)v9djdkxd\ky 


(•) 


denoting 


(10) 
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Figure  1:  Example  of  the  sampling  in  the  first  quadrant  of  the  (kx,  ky)  plane.  Only  the  first 
12  values  of  horizontal  wave  number  k  are  shown,  corresponding  to  m  =  0, 1, 2, ... ,  11.  The 
logarithmic  spacing  in  k  is  already  evident.  The  gridding  is  uniform  in  kx.  Only  12  modes 
need  be  computed  to  support  this  constellation  of  sample  points  in  ( kx ,  ky)  space. 


This  assumes  a  two-sided  spectrum  in  kx.  The  Jacobian 

-a 

arises  from  the  conversion  of  dujdO  to  dkxd\ky\. 


(11) 


Due  to  horizontal  isotropy,  the  computational  effort  of  calculating  the  modes  can  be 
significantly  reduced  over  that  involved  in  a  standard  sampling  of  the  ( kx .  ky)  plane  by  only 

computing  the  modes  on  a  grid  of  k  =  \Jkx  +  ky  values.  The  sampling  here  is  also  chosen 
to  have  logarithmic  spacing  with  the  first  increment  in  each  decade  equal  to  about  10% 
over  the  first  value.  This  puts  about  25  sample  points  in  each  decade.  Thus,  the  horizontal 
wave  number  spacing  is  given  as  km  =  korm,  m  =  0, 1,2, . . .  with  r  =  100  01.  The  lowest 
mode  wave  number  is  set  as  ko  =  2n /L^ho-,  where  L^ho  is  the  length  scale  of  this  mode. 


The  conversion  from  sampling  in  ( kx ,  ky )  to  sampling  in  (k,kx)  on  annuli  of  constant 
k  results  in  non-uniform  contributions  to  the  simulation  along  the  ky  axis  (Fig.  1).  Here, 
the  sample  ( kx,ky )  locations  for  the  first  12  “rings”  (m  =  0, 1,  2, . . . ,  11)  are  shown,  and 
Lkho  =  1  x  105  nr  and  A  kx  =  2n /Lx  with  Lx  =  lx  106  m.  Only  the  (kx,  ky)  sample  locations 
corresponding  to  positive  k‘y  =  k 2  —  kx  are  valid.  This  results  in  sparse  samping  near  the  kx 
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Figure  2:  Examples  of  the  differential  areal  tiling  based  on  the  (k,  kx)  sampling  scheme. 
The  example  presents  only  the  sampling  in  the  first  quadrant  of  the  (kx,  ky)  plane. 


axis.  This  scheme  is  designed  for  accurate  simulations  along  the  x  axis,  but  not  the  y  axis. 

The  differential  variance  contributions  defined  by  Eq.  9  result  in  a  non-uniform  kx  — 
ky  tiling  of  the  ( kx,ky )  plane  (Fig.  2).  A  small  correction  to  the  differential  variance 
contribution  (Eq.  9)  is  required  for  the  first  ring,  or  differential  elements  adjacent  to  the  kx 
axis.  If  m  =  0  or  ki  <  kx,  then  d\ky\  =  yjkj  —  kx.  This  essentially  models  the  contribution 
as  coming  from  a  triangle  with  an  apex  at  the  origin  (instead  of  a  rectangle). 

The  group  speed  dio/dk  is  obtained  directly  using  the  Hellmann-Feynman  theorem  [2,3]. 
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Using  Eq.  3, 


div  lk 
dk  2  uj  X 


2  _  -  — 
X  dk 


(12) 


The  term  dX/dk  can  be  computed  from  Eq.  2.  Differentiating  both  sides  of  Eq.  2  by  A;  yields 

(13) 


2k  dX  , 

^-p'p=Tk'p 


Multiplying  both  sides  by  ijj  and  integrating  over  depth  yields 

2k  j  ip2(z)  dz=^J  ( n 2  -  f2)p>2(z)  dz 
using  the  definition  of  ip.  However,  from  Eq.  4,  f  ( n 2  —  f2)(p2  dz  =  1.  Hence 


(14) 


dco 

dk 


k_ 
uX 
u2-f2 
u>k 


1-y  f<p2dz 


l-{u2  -  f2)  /  ip2  dz 


(15) 

(16) 


Thus,  once  the  modes  and  the  eigenfrequencies  uj  have  been  computed,  the  group  speed 
duj/dk  can  be  determined  directly.  No  WKB  approximations  are  required  in  this  step. 

To  this  point,  any  spectrum  could  be  used.  We  follow  H97  who  utilize  a  GM  spectrum 
modified  so  as  not  to  include  WKB  approximations.  This  model  decomposes  into  frequency 
and  vertical  mode  number.  In  this  representation  the  spectrum  is  separable.  H97  defines 
the  spectrum  as 

S(uj,j)  =  (BEGM)(n0B)2H(j)-^-Vuj2  -  /2  (17) 


7 r  ur 


where 


and 


tiqB  =  J  n(z)dz 


(18) 

(19) 


The  constant  Hq  normalizes  to  unity  the  summation  of  H(j).  The  WKB  limitation  is 
removed  from  the  traditional  GM  spectral  model.  Density  does  not  enter  this  formulation; 
the  spectrum  is  in  terms  of  energy  per  unit  mass.  Note  BEqm  is  a  single  parameter  and 
j*  is  the  traditional  characteristic  (vertical)  mode  number.  Both  BEqm  and  j*  must  be 
specified  in  each  simulation.  GM  neglected  the  contribution  of  the  deeper  part  of  the  ocean 
in  the  integral  of  Eq.  18,  resulting  in  a  significant  underestimate  of  n^B.  If  a  corrected  value 
of  tiqB  is  used,  both  BEqm  and  j*  shoidd  be  appropriately  adjusted.  Note  also  that  Eq.  1 
is  written  as  a  function  of  y.  At  this  time,  we  have  only  implemented  the  code  to  evaluate 
the  field  at  a  single  value  of  yo .  For  convenience,  we  have  selected  the  plane  yo  =  0. 
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buoyancy  frequency  [cph] 


Figure  3:  Buoyancy  profile  used  in  the  tests.  The  profile  is  the  “Munk  canonical  profile” 
modified  by  the  Coriolis  frequency  at  25°. 

3  Tests 


We  present  two  tests  to  demonstrate  both  the  validity  and  viability  of  the  method.  First 
we  examine  the  variance  of  the  simulated  internal  wave  variability.  The  second  explores  the 
distribution  of  energy  in  horizontal  wave  number.  In  both  cases,  the  numerical  results  are 
compared  with  predictions  from  WKB  theory  that  are  readily  calculable. 

In  the  following  tests,  a  test  buoyancy  profile  was  constructed  using  n2(z )  =  + 

/2,  where  raMnnv(z)  =  n§erz!B ,  and  no  =  3  cph  and  B  is  a  vertical  scale  length  here  assigned 
1300  nr.  The  Coriolis  frequency  was  based  on  a  latitude  of  25°.  This  profile  is  shown  in 
Fig.  3. 

With  the  buoyancy  profile  specified,  the  other  parameters  found  in  the  spectrum  (Eq.  17) 
must  be  set.  The  computed  values  of  NqB  and  BEqm  are  6.67  rns^1  and  8.32  x  10-2  vice 
the  canonical  values  of  6.81  nrs-1  and  8.20  x  10~2.  The  differences  are  due  to  numerical 
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Figure  4:  Color  contour  plot  for  a  typical  single  realization  of  internal  wave  displacements. 
A  range  of  100  km  is  plotted  (corresponding  to  an  n/ f  aspect  ratio  of  20).  iseed  was  — 1 
and  t  was  0.0. 

error  in  the  quadrature  of  n(z).  The  low-mode  cutoff  was  assigned  the  value  of  the  canonical 
model,  j*  =  3.  The  computational  domain  was  1000  km  in  range  and  5.5  km  in  depth  with 
sampling  of  8192  points  in  the  horizontal  and  512  points  in  the  vertical.  In  the  program 
a  base  2  FFT  is  used  in  the  horizontal.  This  can  be  changed.  A  vertical  sampling  of  512 
steps  was  chosen  for  convenience.  This  may  be  increased  for  larger  j *.  The  total  number 
of  internal  wave  vertical  modes  was  chosen  to  be  80. 

A  portion  of  a  typical  realization  is  shown  in  Fig.  4.  Note  that  because  of  the  aspect 
ratio  the  contours  are  distorted.  In  reality  the  displacement  contours  form  pancake-like 
structures  stretched  out  in  the  horizontal. 

VARIANCE 

The  simplest  test  is  to  estimate  the  displacement  variance  for  a  given  depth  and  compare 
with  WKB  theory.  The  variance  given  by  Munk  [6]  is 

(C2)wkb  =  \ BEgm (20) 
n(z) 

Multiplying  both  sides  by  n(z)  removes  the  depth  dependence  from  the  Munk  variance,  and 
one  obtains 

n(~)(C2)wKB  =  \ BEgmit-qB  (21) 

which  we  compare  to  n(z)((2).  The  comparisons  between  Eq.  21  and  the  result  from  five 
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Figure  5:  Five  simulation  realizations  of  n(£2)  using  iseed  =  — 1,  — 2,  — 3,  — 4,  — 5.  The 
vertical  line  is  the  WKB  result  at  0.277.  This  figure  tests  Eq.  21. 

simulations  are  shown  in  Fig.  5.  The  simulations  all  used  time  t  =  0,  and  random  number 
generator  seeds,  iseed,  of  —1,  —2,  —3,  —4,  and  —5.  Over  mid-depth,  the  simulated  results 
are  within  ±  20%  of  the  WKB  prediction.  The  comparison  above  2500  m  is  especially  good. 
Very  near  the  surface  and  toward  the  bottom,  WKB  is  not  expected  to  be  valid. 

TOWED  DISPLACEMENT  SPECTRA 

We  turn  now  to  compare  displacement  spectra  of  (one-dimensional)  horizontal  wave 
number.  Referred  to  as  “towed,”  these  spectra  may  be  estimated  from  measurements  using 
horizontal  tows  of  sensors  that  are  suspended  at  fixed  depths  from  a  ship  steaming  in  one 
direction.  The  displacement  may  be  inferred  from  temperature  measurements,  for  example, 
by  using  the  vertical  temperature  gradient  and  assuming  a  linear  relationship. 

From  the  simulations,  we  have  the  displacement  data  directly.  A  sample  spectrum  is 
shown  in  Fig.  6  from  one  of  the  simulations.  There  is  considerable  variability  at  low  wave 
number.  The  theoretical  spectrum  can  be  obtained  from  Eq.  1  and  is 

TSc(kx)  =  S^jtfV9T%{z)  d\ky\  (22) 

j 

Both  spectra  generally  follow  a  k~2  dependence  at  high  wave  number.  This  is  expected 
from  WKB  theory  where  for  convenience,  we  obtain  predictions  by  evaluating  numerically 
an  expression  (Eq.  19)  from  Levine  et  al.  [5].  This  expression  is  complicated  but  may  be 
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Figure  6:  Comparisons  of  horizontal  “towed”  spectra  at  depth  700  m.  A  single  realization  is 
used  (iseed  =  —1).  The  sample  variance  was  92.74.  The  autospectrunr  was  estimated  using 
the  multitaper  method  with  8  tapers  and  a  tinre-bandwidth  product  of  4.  This  is  shown  in 
blue.  The  estimated  integral  of  the  sample  autospectrunr  was  93.77.  The  full  theoretical 
expression  (using  101  modes)  is  shown  in  red  —  this  is  Eq.  22.  The  discretization  of  the 
sampling  in  wave  number  space  is  evident  as  a  “stair-step”  appearance  in  the  theoretical 
curve  at  higher  wave  numbers.  The  predicted  variance  based  on  an  estimated  integral  of 
this  curve  is  83.97.  Additionally,  the  Levine  et  al.  [5]  approximate  expression  is  shown  with 
a  solid  black  line.  This  figure  tests  Eq.  22  and  Eq.  23. 


interpreted  by  considering  the  approximate  formula  given  by  Desaubies  [1]  for  high  wave 
number: 


TSc(fa)  =  *  (4)3  BEaJ-d-  (in  ’j  -  k-2  (23) 

There  is  an  n(z) -1  dependence  in  depth  and  a  k~ 2  dependence  in  wave  number. 

We  expect  from  Fig.  5  that  spectra  near  the  surface  and  bottom  will  not  follow  the 
WKB  behavior.  We  select  a  mid-depth  of  700  nr  for  the  comparison.  The  sample  spectrum  is 
shown  in  Fig.  6  along  with  the  WKB  prediction.  The  model  and  data  spectra  are  consistent 
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with  each  other  and  confirm  the  horizontal  wave  number  behavior  in  the  simulation. 

Both  the  total  variance  and  the  horizontal  spectral  behavior  of  the  simulations  are 
consistent  with  the  underlying  model.  Each  simulation  is  an  independent  realization.  The 
variability  observed  between  simulations  is  typical  for  Monte  Carlo  experiments  and  demon¬ 
strates  the  acceptability  of  the  method. 


4  Summary 


This  report  describes  a  numerical  simulator  for  creating  two-dimensional  range/depth  slices 
of  sound  speed  variability  associated  with  ocean  internal  waves.  The  simulator  is  both 
efficient  and  realistic.  Designed  for  studies  of  internal  wave  effects  on  nreganreter  acoustic 
propagation,  it  may  also  be  used  for  short  and  intermediate  ranges.  The  algorithm  does 
not  depend  on  WKB  approximations,  and  therefore  remains  accurate  in  shallow  water  or 
near  the  ocean  surface.  Combined  with  the  appropriate  deterministic  background  profiles, 
the  statistical  model  for  the  internal  waves  produces  realizations  that  are  true  to  ocean 
measurements  while  the  method  is  easy  to  implement  with  low  computational  cost. 

The  method  also  can  be  extended  to  simulations  wherein  an  accurate  random  wave  field 
is  needed  for  arbitrary  values  of  x  and  y  —  i.e.,  fully  3-D  simulations,  not  only  simulations 
along  a  section  y  =  j/o-  In  such  cases,  the  tiling  in  the  ( kx ,  ky)  domain  should  be  square,  so 
that  dkx  and  dky  are  the  same  everywhere.  Gridding  of  the  kx  and  ky  axes  then  becomes 
uniform.  Efficiency  is  retained  as  the  contributions  {k,  dwjk/dk,  (pjk(z)}  are  assigned  to  be 
the  contributions  from  the  annulus  for  k!  closest  to  the  point  (kx,ky). 
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Appendix:  Pseudo-code 

The  pseudo-code  for  an  implementation  of  the  method  (Fig.  7)  and  input  parameters  (Ta¬ 
ble  1)  are  given.  Where  appropriate,  these  parameters  were  used  in  the  results  described 
in  this  report. 

loop  over  kh  {  kh  =  kO*r**m;  m=0,  1,  . . . ,  nkh  } 

compute  modes,  frequencies,  and  group  velocities  for  each  kh 

loop  over  kx  {  kx  =  0,  1/Lx,  ...,  Nx/(2Lx) 
zero  accumulator  at  kx  and  depth  z 
compute  d|ky| 

loop  over  j  {1 ,2, . . . JMAX} 

compute  variance  Var  =  S (omega, j)*vgp/(pi*kh)*dj*dkx*d|ky | 
generate  G+  from  CN(0,  Var) 
generate  G-  from  CN(0,  Var) 

loop  over  z 

accumulate  [G+*e~ (i*omega*t)  +  G-*e~ (-i*omega*t)] *mode_j (z) 
endloop  (z) 
endloop  (j) 
endloop (  kx  ) 
endloop (  kh  ) 

loop  over  depth  z 

zeta(z)  =  IDFTf  accumulant(kxjz)  } 
endloop  (z) 

Figure  7:  Pseudo-code  for  a  simulator. 
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Parameter 

Description 

Current  Value 

Pseudo-code 

Lx 

x-dimension  (total  range) 

1  X  106  nr 

Lx 

Nx 

number  of  points  in  x-direction 

8192 

Nx 

Nz 

number  of  points  in  z-direction 

512 

3  max 

number  of  vertical  displacement  modes  j 

80 

JMAX 

t 

evolution  time  for  the  realization 

time  [s] 

t 

Nk 

number  of  horizontal  wavenumber  modes  m 

80 

nkh 

LkhO 

length  of  largest  horizontal  mode 

1  X  105  nr 

27r/k0 

iseed 

initializes  random  number  generator 

integer  (<  0) 

3* 

vertical  displacement  mode  bandwidth 

3 

in  SO 

f 

inertial  frequency,  47 rD  sin(latitude)  [rad/s] 

latitude  dependent 

in  SO 

n 

rotation  rate  of  Earth’s  surface 

1  cycle/day 

BEgm 

strength  [nr] 

8.32  x  10“2 

in  SO 

n(z ) 

buoyancy  frequency  profile  [rad/s] 

profile  dependent 

in  SO 

r 

log-step  increment,  horizontal  wavenumbers 

100.04 

r 

Table  1:  Parameters  required  in  the  algorithm. 
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